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A  GENERAL  METHOD  AND  FORTRAN  PROGRAM 
FOR  THE  DESIGN  OF  RECURSIVE  DIGITAL  FILTERS 


INTRODUCTION 

A  very  tedious  problem  often  confronting  engineers  is  the  design  of  a 
frequency  selective  digital  filter  that  meets  a  desired  set  of  specifications.  The 
task  is  especially  laborious  and  the  theory  quite  complex  when  the  filter  has  to 
be  designed  to  satisfy  very  demanding  requirements.  Considerations  such  as 
these  have  often  caused  the  designer  to  settle  for  a  filter  with  characteristics 
further  from  the  ideal  rectangular  type  response  than  is  desired.  Described 
herein  is  a  general  method  for  the  design  of  digital  filters  and  a  versatile 
FORTRAN  program  (listed  in  the  appendix)  that  will  carry  out  the  design  of  any 
low-pass,  high-pass,  bandpass,  or  band-reject  digital  filter  in  any  of  the 
standard  forms  (i.e. ,  Butterworth,  Chebyshev,  or  elliptic).  The  user  need  only 
provide  the  program  with  a  set  of  structured  specifications,  and  as  output  he 
obtains  the  transfer  function  of  the  minimum  order  filter  that  meets  his  require¬ 
ments.  For  convenience,  the  transfer  function  is  expressed  in  a  form  that 
allows  immediate  implementation  is  a  cascade  type  realization.  The  user  also 
has  the  option  of  obtaining  plots  of  the  frequency  response  (both  magnitude  and 
phase)  and  the  unit  sample  response  of  the  filter. 

Complicated  design  nroblems  will  no  longer  demand  so  much  time  and 
effort,  as  has  especially  been  the  case  with  elliptic  filters.  Moreover,  it  is 
suspected  that  elliptic  filters  may  be  preferred  now  for  many  applications  since 
they  normally  require  a  lower  order  filter  than  that  required  by  the  currently 
more  common  Butterworth  and  Chebyshev  filters.  The  advantage  lies  in  the 
fact  that  a  lower  order  implies  fewer  computations  in  the  recursive  scheme  used 
in  implementing  the  filter.  The  principal  argument  against  elliptic  filters  in  the 
past  has  been  the  time  consuming  and  complicated  design  process. 


In  this  report  instructions  for  the  use  of  the  computer  program  are  dis¬ 
cussed  first  aad  then  its  use  by  example  is  illustrated.  Discussion  of  the 
underlying  theory  is  deferred  to  the  latter  part  of  the  report  since  its  compre¬ 
hension  is  not  recessary  in  order  to  successfully  use  the  program  and  the 
resulting  filter. 
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USE  OF  THE  COMPUTER  PROGRAM 

One  data  card,  which  provides  the  program  with  the  desired  set  of  speci¬ 
fications,  should  be  used  for  each  filter  to  be  designed.  The  structure  of  the 
data  card  is  shown  in  table  1.  Also,  figure  1,  which  illustrates  how  the  various 
input  parameters  should  be  interpreted  on  plots  of  typical  frequency  curves, 
should  be  helpful  to  the  user.  More  detailed  descriptions  of  the  input  parame¬ 
ters  and  other  Special  notes  now  follow. 


Table  1.  Data  Card  Structure 


Input  Variable 

Format 

Columns 

Brief  Description 

SRATE 

F8.0 

1-8 

Sampling  rate 

FCLOW 

F8.0 

9-16 

Lower  cutoff  frequency 

FCHIGH 

F8.0 

17-24 

Higher  cutoff  frequency 

RIPPLE 

F6.0 

25-30 

Passband  ripple 

FSLOW 

F8.0 

31-38 

Lower  stopband  boundary 

FSHIGH 

F8.0 

39-46 

Higher  stopband  boundary 

STPLVL 

F4.0 

47-50 

Minimum  stopband  attenuation 

IKIND 

11 

51 

Kind  of  i'ilter  desired  (i.e., 
Chebyshev,  Butterworth,  elliptic) 

ITVPE 

11 

52 

Type  of  filter  desired  (e.g. , 
low-pass,  or  high-pass.) 

NPOLES 

12 

53-54 

Order  desired;  set  equal  to  /,ero 
for  calculation  of  order 

IP  LOT 

11 

56 

Set*0  for  plot  of  magnitude 
of  frequency  response 

NPTS 

14 

57-60 

Number  of  points  in  unit  sample 
response 

FLOW 

FI  0.0 

61-70 

Lowest  frequency  in  plots 

FHIGH 

FI  0.0 

71-80 

Highest  frequency  in  plots 

2 
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The  first  Input  parameter,  SR  ATE,  is  the  sampling  rate  and  Is  normally 
'  the  number  of  data  samples  per  second.  However,  any  time  unit  can  be  used, 
and  in  fact,  one  might  choose  a  fictitious  value  here  with  the  intention  of  con¬ 
structing  a  filter  whose  cutoff  is  a  certain  percentage  of  the  Nyqulst  rate.  For 
example,  a  low-pass  filter  with  a  cutoff  that  is  one-quarter  of  Nyqulst  could  be 
obtained  by  choosing  SRATE  to  be  200  Hz  and  the  cutoff  to  be  25  Hz.  Note  that 
this  same  filter,  if  applied  to  data  sampled  at  some  other  rate,  would  no  longer 
have  cutoff  at  25  Hz,  but  its  cutoff  would  shift  to  one-quarter  of  the  new  Nyqulst 
rate. 


The  quantities  FCLOW  aud  FCHIGH  are  the  cutoff  frequencies  of  the  filter,  | 

and  must  be  given  in  the  same  units  as  SRATE.  Note  that  if  a  low-pass  or  bigh- 
pass  filter  is  desired,  there  is  only  one  cutoff  frequency  and  this  value  should 
be  read  into  FCLOW. 

I 

* 

The  quantity  RIPPLE  is  the  maximum  passband  ripple  (in  decibels!  to  be 
allowed.  In  the  case  of  a  Butte rworth  filter,  where  the  response  is  monotonic, 

RIPPLE  is  not  used.  Note  also  that  RIPPLE  shoula  always  be  a  positive 

quantity.  '  j 


The  parameters  FSLOW  and  FSHIGH  are  similar  to  FCLOW  and  FCHIGH, 
but  they  mark  the  boundaries  between  stopband  and  transition  band.  Here  too 
the  quantity  FSHIGH  is  not  used  in  the  low -pass  or  high-pass  case. 


The  next  specification,  STPLVL,  is  the  minimum  stopband  attenuation 
(in  decibels)  the  user  will  tolerate.  This  should  also  always  be  a  positive 
quantity. 

i 

| 

The  value  1,  2,  or  3  of  the  variaule  IK3ND  determines  which  kind  of 
filter  will  be  designed,  that  is,  Chebyahev,  Butterworth,  or  elliptic,  respec¬ 
tively. 


Similarly,  the  value  1,  2,  3,  or  4  of  ITYPE  determines  the  shape  of  the 
filter  that  will  be  designed,  that  is,  low-pass,  high-pass,  bandpass,  or  band- 
reject,  respectively. 
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The  value  of  NPOLES  should  normally  be  read  in  as  zero,  and  the  correct 
order  will  be  calculated.  However,  the  user  may  (except  in  the  elliptic  case) 
choose  a  positive  integer  here  if  he  desires  a  filter  of  some  specific  order. 
This  integer  should  be  the  order  of  the  basic  low -pass  structure  desired  before 
transformation.  In  converting  a  lew -pass  filter  to  either  a  bandpass  or  band- 
reject  filter,  the  method  doubles  the  number  of  poles;  thus,  to  get  a  band-reject 
filter  of  order  18,  NPOLES  should  be  set  equal  to  9.  If  NPOLES  is  given  as 
nonzero,  the  quantities  FSLOW,  FSHIGH,  and  STOLVL  are  not  used  since  the 
given  order  determines  these  quantities  uniquely.  NPOLES  is  restricted  to  be 
less  than  or  equal  to  20,  since  it  is  anticipated  that  filters  of  these  orders  can 
meet  any  reasonable  set  of  requirements.  However,  if  necessary,  one  could 
merely  increase  the  array  dimensions  used  in  the  program  in  order  to  obtain 
filters  of  higher  order. 


The  remaining  input  parameters  control  the  plotc  produced  by  the  program. 
If  tii8  user  does  not  have  access  to  the  Stromberg  Carlson  4060  Integrated 
Graphics  System,  he  should  remove  the  corresponding  code  or  replace  it  by  a 
code  that  is  compatible  with  his  plotter. 


If  for  some  reason  the  user  doe?  not  want  any  of  the  three  available  plots, 
he  should  set  the  parameter  IPLOT  to  zero.  The  number  of  points  desired  in 
the  uv  it  sample  response  should  be  read  into  NPTS.  This  quantity  must  not  be 
greater  than  1000  and  has  a  default  value  of  100  if  the  corresponding  field  in  the 
data  card  is  left  blank.  The  values  given  to  the  quantities  FLOW  and  FHIGH 
merely  restrict  the  plots  of  magnitude  and  phase  to  the  range  between  these  two 
frequencies.  411  information  can  be  seen  by  plotting  from  zero  to  Nyquist,  since 
the  curve  is  periodic  with  period  SKATE  and  is  symmetric  about  Nyquist.  How¬ 
ever,  one  can  expand  a  small  frequency  range  for  closer  examination. 


If  more  than  one  filter  is  desired,  a  similar  data  card  for  each  should  be 
used  and  a  blank  card  should  be  placed  last. 


The  following  typical  example  Illustrates  the  points  just  discussed.  Sup¬ 
pose  that  an  elliptic  bandpass  digital  filter  with  cutoffs  at  2000  Hz  and  3000  Hz 
is  desired.  Assume  that  the  sampling  rate  is  10,000  Hz,  that  the  filter  is  to  be 
at  least  30  dB  down  by  1800  Hz  and  3200  Hz,  and  that  the  tolerable  ripple  in  the 
passband  is  0.5  dB.  The  fo^owing  is  a  list  of  the  proper  input  parameters  for 
tnis  problem: 
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SRATE 

10000.0 

FCLOW 

2000.0 

FcraGH 

3000.0 

RIPPLE 

0.5 

FSLOW 

1800.0 

FSHIGH 

3200.0 

STPLVL 

30.0 

IKIND 

3 

ITYPE 

3 

NPOLES 

0 

IP  LOT 

1 

NPTS 

200 

FLOW 

0.0 

FfflGH 

5000. 0 

Figures  2  through  5  show  the  printed  and  plotted  outputs  for  the  example 
In  question.  Note  that  the  range  of  the  magnitude  plot  is  from  0  dB  to  -50  dB, 
which  can  easily  be  increased  if  desired. 
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InPUl  SPECIFICATIONS 

SKATE  s  lOOOO.OO 

FCLpW  :  2900,00 

FCrtlGH  :  3000,00 

RIPPLE  =  .500 

FSU/»  s  1800.00 

FSHJGH  =  3200,00 

STPLVL  =  30,00 
lftlNO  =  3 

lTTPt  =  3 

NPOLES  =  q 

IPLOT  =  1 

InE  MINIMUM  uROErt  FIlTEm  NHICI*  MEETS  THESE  SPECIFICATIONS  IS  ORDER  8 


InE  IrtANsFErt  FUNCTION, H(2) ,  oF  THE  OtSIREO  FILTER  CAN  BE  EXPRESSED 
in  TrtE  FuLLOnINb  FORM 


M(Z> 


t-NsT  « 


2«*2»8(1)»ZHI<2> 

2*«2*C(1)*2+l<2I 


2«»2*B(3)«2*B(»t) 

2»»?>C(3)»Z*C(4) 


•hEKc  Tut  COlFFUIEUS  «,<C  GjVEH  8Y 


CNST  =  .313160351429645-001 

Cl 

11  = 

—  29j707o452a76964000 

6( 

11  = 

. 1406719642914184001 

Cl 

21  = 

. 744807 7036l740=4QQo 

d( 

21  = 

.1009000030000004001 

C  { 

3)  = 

.  24370  7u4528769iI  4000 

d( 

31= 

-. 146=71 954291 4 184001 

Cl 

41  = 

,74h807  70361740=4000 

b( 

41  = 

.1000000000000004001 

Cl 

51  = 

- , 606822*5784317=4000 

8< 

51  = 

.6926035336538994000 

cl 

61  = 

.93041 1h926636oo*OOo 

B( 

61  = 

,1000000000000004001 

c( 

7|  = 

. 60 o822457643l 7=4000 

o( 

71  = 

-.8925035336538994000 

Cl 

81  = 

.9304114926636604000 

B( 

81  = 

,1000000000000004001 

Im£  TKAnsFEH  FUNCTION  IS  EXPhLSSc.0  III  The  ABuVE  form  so  that 

I  HE  LuEFl  ICIlNTS  can  BE  oSEO  AS  »IVEN  IN  A  CaSCAOE  IMPLEMENTATION  OF  THE  FILTER 


Figure  2.  Sample  Printout 


Z**24BtN-l)*Z4B(Ml 

Z**24C(N-1»*Z*CIN) 
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360. 


Frequency  Characteristic  (Phase) 


Figure  5.  Unit  Sample  Response 
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Ti!EOFv  "10GRAMMING  METHOD 

PRELIMINARIES 

Each  of  the  major  steps  presented  in  the  following  synopsis  of  operations 
performed  by  the  program  will  be  discussed  in  greater  detail  in  succeeding 
sections: 

a.  Certain  critical  frequencies  in  the  digital  z-plane  are  given  as  input 
parameters.  These  must  be  transformed  to  continuous  s-plane  critical  fre¬ 
quencies  since  much  of  the  design  process  is  carried  out  in  the  s-plane. 

b.  The  minimum  order  filter  of  the  appropriate  kind  that  will  meet  the 
specifications  is  computed. 

c.  The  s-plane  pole-zero  pattern  of  a  unity  bandwidth  low-pass  analog 
filter  of  proper  order  is  determined. 

d.  These  s-plane  poles  and  zeros  are  mapped  to  the  z-plane  poles  and 
zeros  of  the  required  filter. 

e.  The  poles  and  zeros  of  the  desired  transfer  function  determine  the 
function  up  to  some  constant  factor;  the  required  constant  is  obtained  by  forcing 
the  maximum  value  of  the  amplitude  characteristic  to  be  unity  (0  dB). 

f.  Next,  the  transfer  function  is  manipulated  into  a  form  that  can  be 
immediately  implemented  in  cascade  form. 

g.  The  magnitude  and  phase  of  the  transfer  function  evaluated  on  the  unit 
circle  in  the  z-plane  (i.e.,  the  frequency  response)  is  plotted  in  the  frequency 
range  specified  by  the  user. 

h.  Finally,  the  unit  3ample  response  of  the  desired  filter  is  plotted. 


11 
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A  few  facta  from  Alter  theory  are  now  presented  as  background  for  later 
discussion.  The  three  kinds  of  fillers  considered  in  this  report  have  magnitude 
characteristics  in  the  analog  plane  defined  by  the  following  equations : 


i 


Butterworth 

|HflO)|2  “  “S 

(la) 

i 

1  +fi 

1 

h 

1 

i 

| 

Chebyshev 

|HW|\Jv> 

(lb) 

elliptic 

N«)l2--  22“' 

(lc) 

j 

1  r€  in(Q) 

where 


H(* )  -  the  transfer  function  of  the  unity  bandwidth  low-pass  filter 
(1  =  radian  frequency 

n  =  the  order  of  the  filter 

€2  =  10R/1°  -  1,  where  R  is  the  amount  of  passband  ripple 

(in  decibels) 

V  (•)  =  Chebyshev  polynomial  of  order  n 

\  t 


*.»  - 


sn 


sn 


n  odd 


r  K^i>  -i  i 

Ln  KCkT“  (ft’k)’\r 

r  K<ki>  -i  i 

[WJ  +  n  sn  (fi,  k),  kj ,  n  even, 


where 


sn(* ,  •  )=  Jacobian  elliptic  function 

K(. )  =  complete  elliptic  integral  of  the  first  kind 


i 

r 
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s  s 

1/ft  ,  ft  being  the  start  of  the  stopband 
j  2  S/20 

C/vA  -  1,  with  A  =  10  and  S  equal  to  the  minimum 
stopband  attenuation  (in  decibels). 


Two  types  of  mappings  that  play  a  major  role  in  the  method  will  now  be 
presented.  The  first  is  a  map  that  converts  a  unity  bandwidth  low-pass  analog 
filter,  say  H(s),  into  an  analog  filter  of  another  type  with  different  cutoff(s). 
If  H(s)  is  such  a  filter,  then  v-e  see  that 


and 


H(s/ftC) 

H(ftC/s) 


H 


fs2 


(J 

is  low-pas?  with  cutoff  ft 

Q 

is  high  -pass  with  cutoff  ft 

c  c 

is  bandpass  with  cutoffs  ft.  and  ft2 


(2a) 

(2b) 


(2c) 


H 


s<«2  " 


2  nc 

8  +  fll  n2. 


c  c 

is  band-reject  with  cutoffs  ft  and  ft  .  (2d) 

±  & 


A  second  mapping,  called  the  bilinear  transformation,  maps  the  s-plane 
to  the  z-plane  and  can  be  used  to  convert  an  analog  filter  to  a  digital  filter.  The 
mapping  is  given  by 


Note  that  the  imaginary  axis  in  the  s-plane  is  mapped  to  the  unit  circle  in  the 
z-plane  and  that  the  left  half  of  the  s-plane  is  mapped  inside  the  unit  circle 
in  the  z-plane.  Thus,  a  stable  analog  filter  will  be  mapped  to  a  stable  digital 
filter  and  the  digital  frequency  response  from  zero  to  Nyquist  will  take  on 
exactly  the  same  values  as  the  analog  frequency  response  from  zero  to  infinity. 
Also,  note  that  the  aliasing  problem  inherent  in  filters  designed  by  the  method 
of  impulse  invariance  is  not  present  here  since  the  mapping  is  invertible.  It  is 
often  stated  that  this  method  has  the  drawback  of  warping  the  frequency  scale. 
However,  no  real  problem  exists,  since  the  critical  frequencies  in  a  design 
problem  can  be  "prewarped,  "  as  described  in  the  next  section. 
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In  the  literature,  the  terms  Butte rworth,  Chebyshev,  and  elliptic  are  used  ‘ 
in  the  description  of  low-pass  analog  filters.  In  this  report  these  terms  are 
also  used  in  describing  filters  of  other  types  that  have  been  obtained  by  the  pre¬ 
vious  two  mappings. 


ANALOG  CRITICAL  FREQUENCIES 

The  given  digital  critical  frequencies  must  be  transformed  to  the  analog 
plane  in  such  a  way  that  when  they  are  later  mapped  back  to  the  z-plane  by  the 
bilinear  transformation  they  will  be  mapped  to  the  proper  values.  This 
"prewar  '  ng"  of  the  critical  frequencies  is  carried  out  by  the  map 

fi  =  tanj~^  ,  (4) 

where 


fl  =  continuous  radian  frequency 
u  =  digital  radian  frequency 
T  =  time  between  samples. 

The  following  argument  will  convince  the  reader  of  this  result.  Using  the  in¬ 
verse  of  the  bilinear  transformation  (3),  we  can  relate  any  z-plane  point  to  its 
corresponding  s -plane  point,  i.e., 

z  -  1 

a  =  —  * 

For  z  =  eJwT,  at  which  points  we  obtain  the  digital  frequency  response,  the 
corresponding  s -plane  point  is 


e3wT- 1  „ 

8=T3f~~  =  jtan 

eJ  + 1 


But  the  analog  frequency  response  is  obtained  by  evaluating  the  transfer  function 
at  s  =  jfi .  Thus,  we  have  the  correspondence  dictated  by  (4).  Now  when  we 
apply  the  bilinear  transformation  to  cur  s-plane  filter,  we  can  be  sure  we  will 
end  up  with  a  filter  with  proper  digital  critical  frequencies. 
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DETERMINATION  OF  MINIMUM  ORDER 

We  will  first  solve  the  problem  for  the  low -pass  case  and  then  show  how 
any  of  the  other  cases  can  be  reduced  to  an  equivalent  problem. 

c  c 

Let  oj  and  fi  be  the  desired  digital  and  analog  radian  cutoff  frequen¬ 
cies,  respectively,  and  let  «s  and  (Is  be  the  frequencies  at  which  the  stopband 
begins.  By  (4),  these  frequencies  are  connected  by  the  relations 


f°  and  f3  are  the  given  input  specifications  in  hertz. 

Next,  define  the  "transition  ratio,  "  an  important  factor  in  determining  the 
required  order,  by 


nc 


Using  the  expressions  (la)-(lc)  given  previously  for  the  three  magnitude 
characteristics,  we  can  determine  N,  the  minimum  order  that  will  satisfy  the 
filter  requirements.  Omitting  details  of  the  algebra,  we  now  present  general 
expressions  for  calculation  of  N. 


If  a  Butterworth  filter  is  desired,  N  can  be  computed  as  the  smallest 
integer  that  is  greater  than 

log10<10S/l°-l> 


2  i°g10  o 


15 


TR  4629 


For  a  Chebyshev  filter,  we  compute  the  required  order  by  finding  the 
smallest  value  of  N: 

T  /lOS/1°  -  l\  1/2 

n<  Ur/io-J  ■ 

Finally,  if  an  elliptic  filter  is  desired,  we  determine  N  by  finding  the 
smallest  integer  that  is  greater  than 


k(Ji  -k|)K(l/fl- 


Kfkj)  K  [yl  -  (1/  )  J 

All  the  parameters  used  in  these  expressions  have  been  defined  in  this  or 
previous  sections. 

The  only  additional  problem  arises  in  the  elliptic  case  and  results  from 
the  fact  that  the  expression  above  will  ordinarily  not  be  an  integer  and  must  be 
rounded  upwards.  Since  it  is  required  that 


k(  /l  -  k ^  )  K(J/  flT) 


k[vi  -(1/n  )  J 

the  value  of  kj  must  be  recalculated  after  N  is  determined.  A  convenient 
formula*  for  recomputation  of  kj  is 


2,0'25  r!+£qi(i+1>T 

L  i=i  J 


1  +  2  2^ q 

i=l 


where 


q  =  exp 


!NnK(7l 

K(k 


-k2) 


Note  here  that  k  =  1/  fi  . 


*R.  M.  Fano,  "A  Note  on  the  Solution  of  Certain  Approximation  Problems 
in  Network  Synthesis, "  J.  Franklin  Institute,  vol.  249,  1950,  pp.  189-205. 
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Now  the  problem  of  determination  of  minimum  order  for  the  low-pass  case 
is  solved.  For  any  other  case  we  define  ftT  differently  but  proceed  exactly  as 
in  the  low -pass  case.  For  the  high-pass  case  the  appropriate  definition  is 


nT  =  IL. 

fiS 


For  the  bandpass  case  chcose 

flT  =  min  (ilj ,  ftj)  , 
and  for  the  band-reject  case  choose 


\ni  fl2/ 


where  the  ft  f  are  given  by 


,„S,2 

(fi.)  -  fij 


<»2 


«;> 


i  =  1,2. 


C  c  s  s 

By  ftj  and  &2  we  mean  the  two  cutoff  frequencies,  and  by  ftj  and  0.%  we 

mean  the  two  boundaries  between  stopband  and  transition  band. 


T 

These  rules  for  determination  of  ft  and,  hence,  N,  for  filter  types 
other  than  low-pass,  follow  directly  from  (2a)-(2d),  the  transformations  that 
convert  a  low-pass  filter  to  one  of  another  type.  Note,  however,  that  in  the 
bandpass  and  band-reject  cases,  the  computed  N  will  equal  the  order  of  the 
basic  low -pass  structure  before  transformation  to  the  filter  of  proper  type. 

The  actual  order  of  the  final  filter  will  be  2N. 


In  the  above  we  have  chosen  the  minimum  order  filter  which  meets  the 
user's  requirements.  It  should  be  mentioned.,  however,  that  the  user  may  find 
another  order  more  appealing  in  certain  special  cases.  For  example,  for  band- 
reject  Chebyshev  or  elliptic  filters  that  are  derived  from  even  order  low-pass 
structures,  we  find  that  the  magnitude  characteristic  is  asymptotic  to  -R  dB  in 
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the  passband,  and  is  actually  near  -R  dB  over  a  large  percentage  of  the  frequency 
range.  Thus,  a  filter  derived  from  an  odd  order  low-pass  structure  may  be 
preferred  here,  since  its  frequency  response  is  asymptotic  to  0  dB,  the  ideal 
value  in  the  passband. 


DETERMINATION  OF  S- PLANE 
POLE -ZERO  PATTERN 

A  lengthy  but  straightforward  task  is  the  computation  of  the  pole-zero 
pattern  of  a  continuous  unity  bandwidth  low-pass  filter  of  order  N.  The  tech¬ 
niques  are  adequately  reviewed  by  Gold  and  Rader*  and  will  not  be  repeated 
here.  However,  it  is  not  necessary  to  keep  track  of  all  poles  and  zeros  since 
they  always  occur  in  complex  conjugate  pairs,  and  it  is  not  necessary  to  com¬ 
pute  any  zeros,  except  in  the  elliptic  case,  since  they  are  always  at  infinity  in 
the  s -plane. 


Also,  the  complete  elliptic  integral  of  the  first  kind,  which  is  needed  in 
the  elliptic  design  problem,  is  given  exactly  by 

1 

=  J  7  2  2  2  1l/2  ’ 

0  [(1  -  t  )(1  -  k  t  )J 

but  in  the  program  is  approximated  by  the  formula^ 


K(m)  =  (i 


a^  +  a^m^  +  •  •  •  + 


V'l)  *(b0*bl"lt"'*  b4ml )  to  (^) 


*B.  Gold  and  C.  M.  Rat  r,  Digital  Processing  of  S'gnals,  McGraw-Hill 
Book  Company,  Inc.,  1969,  ph  48-97. 

^M.  Abramowitz  and  I.  A.  Stegnn,  Handbook  of  Mathematical  Functions, 
Dover  Publications,  1964,  pp.  567-607* 
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aQ=  1.38629436112 

aj  =  0. 09666344259 

a  =  0. 03590092383 

a3  =  0.03742563713 

a.  =  0.01451196212 
4 


b0=°.5 

bj  =  0. 12498593597 
b2  =  0.06880248576 
b  =  0. 0322835534b 

o 

b .  =  0.00441787012 
4 


m  =  k 

m  =  1  -  m. 

* 

The  maximum  error  in  this  formula  is  2  x  10  when  0  <  m  <  1,  in  which  range 
m  will  always  be. 

Also,  other  approximations*  are  used  in  the  calculation  of  elliptic  poles  and 
zeros.  Iterative  and/or  series  calculations  can  be  used  to  replace  these  approxi¬ 
mations  but  the  results  obtained  with  them  were  deemed  sufficiently  accurate. 


CALCULATION  OF  POLE-ZERO  PATTERN 
OF  THE  DESIRED  DIGITAL  FILTER 

Once  the  s-plane  poles  and  zeros  of  a  unity  bandwidth  low-pass  continuous 
filter  of  proper  order  and  kind  have  been  determined,  it  is  possible  to  map  these 
to  the  poles  and  zeros  of  the  desired  digital  filter,  which  may  be  of  any  type  with 
any  cutoff(s). 


Although  the  calculations  vary  by  type  of  filter  desired,  the  method  is 
basically  the  same,  and  will  be  shown  for  the  bandpass  case  only.  If  H(s)  is  the 
transfer  function  of  a  unity  bandwidth  low-pass  continuous  filter  of  order  N,  then 


Ha(s)  =  H 


s2  *  "l 


*Gold  and  Rader,  op.  cit. 
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is  the  transfer  function  of  a  continuous  bandpass  filter  of  order  2N  with  cutoffs 
and  g£. 


Nov  applying  the  bilinear  transformation  (3),  we  see  t>at 

2 

FYili 

ESI" 


Hd(z)  =  Ha 


m)  ♦ 


m  &  ■  <0J 


is  the  transfer  function  of  a  bandpass  digital  filter  of  order  2N  with  cutoffs  at 

c  2  .  -1 

(. 

and 

**2 


c  2  .  -1  /_c\ 

IT  tan  (g^ 


Now  if  we  solve  the  equation 


s  = 


(g) 

Q  « -  ai) 


for  z  in  terms  of  s,  we  obtain 


2  2  1/2 
_  (p  -  1)  *  (s  f  -  4pP  * 

sd  -  1  -  p  * 


where 


v*a°al 


d=  a'  a”  . 
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We  can  merely  substitute  the  previously  computed  s-plane  poles  and  zeros  into 
this  expression  to  obtain  the  poles  and  zeros  of  the  desired  digital  bandpass 
transfer  function  Hq(z).  Note  that  each  point  in  the  s-plane  is  mapped  to  two 
z-plane  points,  and  thus  the  number  of  poles  and  zeros  doubles.  Actually,  the 
computations  are  not  carried  out  for  the  zeros  unless  we  are  working  with  an 
elliptic  filter  since  Hd(z)  will  always  have  N  zeros  at  z  =  1  and  N  zeros 
at  z  =  -1.  This  follows  from  the  fact  that  the  s-plane  zeros  are  all  at  infinity. 

When  we  have  calculated  the  z-plane  poles  and  zeros,  we  have  determined 
the  desired  transfer  function  up  to  some  constant  factor.  In  the  next  section  we 
will  determine  that  constant. 


CALCULATION  OF  THE  CONSTANT 
MULTIPLIER 


Thus  far  we  have  determined  the  transfer  function,  HD(z),  of  the  desired 
filter  as 


N 

><n  <*-*•> 

v'-N— — 

n--» 

i=l 


where 


K  =  constant  multiplier  to  be  determined 
N  =  order  of  filter 
z.  =  calculated  zeros 

l 

p.  =  calculated  poles . 

Knowing  the  value  of  Hq(z)  for  any  z  (z  t  z^  or  p^)  will  enable  us  to  deter¬ 
mine  K  uniquely. 
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and  solve  for  the  quantity  K. 

From  the  theory  of  the  previous  section  it  follows  that  the  proper  values 
for  zc  are 

low -pass  zq  =  1 

high-pass  zq  =  -1 

bandpass  z  =  l — £  ±  j  ,  where  P  =  ft?  »  as  before, 
Ol-p  1+P  16 

» 

band-reject  zq  =  ±  1 . 

For  example,  in  the  bandpass  case  we  know  that  any  point  s0  in  the  s-plane 
will  be  mapped  to  the  point 

2  2  1/2 
(P  -  1)  ±  (sQ  d  -  4p) 

jr  —  —  i  iwi  — — i  ji  ■  . .  • 

o  s  d  -  1  -  p 
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Hence,  the  point  s0  -  0  is  mapped  to 

o  1  +  p  J  1  +  p 

Note  that  in  the  bandpass  and  band-reject  cases  we  have  a  choice  of  two  values 
for  zQ,  since  each  is  an  image  of  the  point  s  =  0. 

No  matter  what  type  filter  we  are  designing,  it  follows  that  we  can  deter¬ 
mine  the  quantity  K  by 


N 


i=l 


In  computing  this  expression,  computational  techniques  make  it  unnecessary  to 
work  with  complex  quantities. 

MANIPULATION  INTO  CASCADE  FORM 

By  combining  the  complex  conjugate  poles  and  zeros  of  the  transfer  func¬ 
tion  Hd(z)  into  second  order  factors  with  real  coefficients,  we  can  write  the 
transfer  function  as  a  product  of  terms  of  the  form 

2 

z  +  b  z  +  b 

i  Ct 

~~2  _  ’ 

Z  +Cj  z  +c2 

There  will  also  be  one  factor  of  the  form 

z  +  b 


if  the  order  of  the  filter  is  odd,  corresponding  to  the  real  pole.  However,  when 
one  actually  implements  a  filter,  some  recursive  scheme  is  needed  to  calculate 
the  present  filtered  data  point.  With  the  transfer  function  in  the  above  form  we 
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can  immediately  write  a  set  of  difference  equations  that  describes  the  filler  in 
the  time  domain,  and  this  set  will  correspond  to  a  cascade  realization.  In  the 
general  case  where  the  transfer  function  is 


HW-ai-Kl 

X(z) 


z  +  b„\/  z2  +  h  z  +  b„ 


z  +  b  z  +  b 
N-l  N 


,z  +  c 


V  \z  +  c2  z  +  c  , 


"TF 


'z  tcN-lz+V 


the  following  set  of  difference  equations  equivalently  describes  the  filter: 
yjin)  ~  x(n)  -  e^fn-l) 
y2(n)  =  yj(n)  +bJy1(n-l) 

Jgfa)  =  y2(n)  -  c2y3(n-l)  -  c3yg(n-l) 
y4(n)  =  y3(»)  +  b2y3(n-1)  +  b3y3(n-2) 
y5(n)  =  y4(n)  -  c4y5(n-l)  -  c,.y5(n-2) 
yg(n)  *  y5(n)  +b4yg(n-l)  +bgy5(n-2) 


yNW  =  '  °N-1  V""  •  CN  yN(n'2) 

y(n)  =  K(yN(n)  +  yN(n-I)  +  bN  yN  (n-2»  , 

where  x(n)  and  y(n)  are  the  filter  inputs  and  outputs,  respectively,  and  yj(n), 
yj(n),  ....  y^fn)  are  intermediate  values  determined  by  the  above  equations. 

The  above  filter  could  also  be  described  by  a  single  difference  equation  of 
order  N,  which  would  be  an  entirely  equivalent  time  domain  representation,  if 
all  coefficients  could  be  represented  exactly.  However,  due  to  the  necessary 
quantization  of  coefficients,  this  so-called  "direct  form"  realization  is  not 
recommended  since  it  is  much  more  sensitive  to  coefficient  accuracy  than  the 
corresponding  cascade  implementation  <5). 
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THE  MAGNITUDE  AND  PHASE 
OF  THE  FREQUENCY  RESPONSE 

To  obtain  the  frequency  response  of  a  digital  filter,  we  must  evaluate  the 
transfer  function  on  the  unit  circle  of  the  complex  z-plane.  In  particular, 
|HD(ej‘‘JT)|>  the  magnitude  of  the  frequency  response,  is  usually  of  primary 
interest. 

Thus  far  we  have  Hd(z)  expressed  as  a  product  of  terms  of  the  form 

2  , 

z  +  b,  z  +  fr 
1  2 

2 

Z  +  c  z  +  c 

I  L* 

The  magnitude  and  phase  of  this  term  for  z  --  eiu^  can  easily  be  expressed  as 
a  trigonometric  function  of  wT  and  can  be  simply  evaluated  for  any  u?.  The 
product  of  the  magnitudes  of  all  such  terms  and  the  constant  multiplier  is  then 
the  magnitude  of  the  frequency  response.  Similarly,  the  sum  of  the  phases  of 
all  such  terms  gives  the  phase  of  the  frequency  response.  The  possibility  of  a 
single  term  of  the  form 

z  +  b 
z  +  c 

poses  no  additional  problem.  In  the  program  we  compute  the  magnitude  and 
phase  as  above,  convert  the  magnitude  to  decibels  and  phase  to  degrees,  and 
then  provide  plots  in  the  frequency  range  desired. 

THE  UNIT  SAMPLE  RESPONSE 

The  rate  of  decay  of  the  unit  sample  response  (i.e. ,  the  response  of  the 
filter  to  a  pulse  of  unit  height  at  time  zero)  of  a  digital  filter  is  often  of  interest 
to  the  designer.  For  this  reason  a  plot  is  provided  if  requested. 

Since  the  transfer  function,  Hq(z),  can  be  viewed  as  the  /-transform 
of  h(n),  the  unit  sample  response,  we  can  write 

h(n)  =  Z"1  [hd(z)J  , 
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or 

h<n)  =~  ^  Hd(z)  zD-1  dz, 

C 

where  C  is  any  contour  In  the  complex  plane  enclosing  all  poles  of  the  integrand. 
Replacing  Hd(z)  by  its  pole-zero  factorization,  we  obtain 


Since  the  pj,  i  =  1,2 . N  will  always  be  simple  poles,  the  residue  at  a 

specific  pj  of  the  integrand  is 


For  the  special  case  n  =  0  there  is  also  a  pole  at  z  =  0,  but  the  complications 
introduced  by  the  additional  pole  need  not  be  considered  since  we  know  by  exam¬ 
ining  (5)  that  h(0)  must  equal  K.  Defining 

Hu  -  v 

i=l  J 

Pr^ - • 

i=l 

i+j 


we  can  then  easily  compute  the  unit  sample  response  as 
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N 

2  P^"1  •  ^ 

i=l 

K  ,  n=G 

Although  the  unit  simple  response  is  a  discrete  function  of  time,  the  plot  pro¬ 
duced  by  the  program  is  a  continuous  function  for  ease  in  viewing. 


SUMMARY 

Butterworth,  Chebyshev,  and  elliptic  digital  filters  are  realizable  networks 
that  approximate  physically  unrealizable  filters  having  rectangular  magnitude 
characteristics.  Each  of  these  filters  come  arbitrarily  close  to  the  ideal  re¬ 
sponse  as  the  order  of  the  approximating  network  becomes  large.  Thus,  the 
question  of  the  order  necessary  to  meet  a  set  of  specifications  is  usually  criti¬ 
cal.  The  computer  program  described  in  this  report  computes  the  minimum 
order  filter  that  meets  the  particular  requirements  of  the  user,  and  then  pro¬ 
ceeds  to  design  that  filter.  The  program  also  produces  plots  of  the  magnitude 
and  phase  of  the  frequency  response  and  the  unit  sample  response  of  the  filter. 
The  reader  need  not  fully  understand  the  theory  presented  in  this  report  in  order 
to  use  the  program.  Finally,  for  those  unfamiliar  with  digital  signal  processing 
techniques,  the  instruction  on  how  the  filter  transfer  function  can  be  implemented 
as  a  set  of  difference  equations  should  be  useful. 
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c 

t ET «tEN  STOPBAND  AND  TRANSITION  BAND 

c 

1 

(IN  CTCLLS/TIHC  UNIT) 

V* 

V. 

r 

alPkVL 

RtAL 

n7-5o 

THE  MINIMAL  STOPBAnD  ATTENUATION! IN  DB) 

s» 

V 

iAINu 

integer 

SI 

THE  MNO  01  FILTER  DESIRED  l.S. 

t 

I  CHCBTsHtT 

c 

2  BUTTErVORTH 

c 

3  ELLIPTIC 

c 

A  TYPE 

INTEGER 

32 

TYPE  OF  FILTER  DESIRED  I.E. 

c 

1  LOV-PASS 

c 

2  HIGH-PASS 

C 

3  BANDPASS 

c 

4  band-reject 

w 

c 

NPOLES 

Integer 

a  3-54 

SETSO  FO«  CALCULATION  OF  REGUlRED  ORDER 

c 

otherwise  setsoroer  oesireo  for  basic 

C 

LOV-F'ASS  STRUCTURE  BEFORE  TRANSFORMATION 

u 

C 

iplot 

integer 

36 

S£T=l  FOR  PLOT  OF  MAGNlTUOE  OF  H(2> 

k  SETSO  FOR  NO  PLOT 


I 


3 

* 


i 

l 


30 


v» 

I'.TEGEP 

o7-6q 

NO.  OF  PTS  In  UNIT  SAMPLE  RESPONSE 
(DEFAULT  VALUE  =  100) 

c 

L 

KlAL 

Ol-7o 

l0«  E‘*P  FREUuEnCT  OF  PLOT 

c 

l 

w 

Hil«h 

Real 

n-Bo 

HIGH  tnO  FREGUENCY  OF  PLOT 

c 

c 

(WTfc— 

—  IF  NpOLES.NL.O  (I.E. 

THE  ORutR  IS  SPECIFIED).  FSLOM.FSHlbH  ANO 

c 

sTPL*L  are 

Not  OSLO 

(UNLESS  ELLIPTIC)  SINCE  THE  GIVEN  ORDER 

Y  uETEuPInES  InEsC  -UA.VTITIES 

Y 

y  IF  YvA-PASS  OR  HlGH-PASS,  FOiIGH  AND  FSHIiH  ARE  NOT  USEO 

Y 

Y  IF  boTTERVORfH  IS  DEslREll.  RIPPLE  IS  NOT  USED 
L 

y  I  HE  .UAJTITIlS  filPPLc.  Ali(J  STPLVL  SHOULD  ALaAYS  BE  GIVEN  .ST.  0 

y 

Y  KF.PEaT  DATA  YARD  OF  IMls  TTPl  FOk  each  filter  to  f>e  DESIbMED 

y  plACy  A  RLAiw,  data  CARO  yAST 

>  ••«**•••«»•••♦••********* 

Y 

Y 


OOubLE  PREy  I  SION  CoEF  <  21) ,  2 ) .  POLES  (10,2)#  POLE2  (40,2).  ZEROS!  10,2) 
UOoplE  PRECISION  y„EFN(F0,2),CNST,y£R0yU0,2) 

OOohlE  PRECISION  EPS2.pI,«CUT1,),COT2,EaPNT,AN6LE,OaHSLE,HSTOP1, 

. »alOP<!,  OHEuA.OHEGy.1,  a.  H»C,U«e.F.S,pRoO,t)IFF»RHAG,  HAGGLE,  THETA. 
2Si<l,Sll0,CNl,CN3,0i-.-r)N»,SN,G 
DIMENSION  a(1000).y(IOoO),PHI(1000),Z(200) 

DIMENSION  .121) 

DIMENSION  m1(1000,,h2(1000) 

COMPLEX  ZPuLESI “0 ) , ZZErOS (40) #£ TEMP (40), BETA (40) 

ComPLEa  PTtMPI.RTtMP2,zIEMH2 

ELlRTC ( T )=l. 38629a Jol 1*OOf  !•  ( .096bb344,:S900*T» ( , 03SR009238300F 
iT»i.'337A2bo37Ijn0»l».0lA5U962l200)  ))F0LOu(I,0/T)»(  ,S*T» 
2(.l2A9dS93a97U0*T«(.06Ho02»857600FT«(,0332S3Sb3AbD0FT*.004A178Y012 
300)))) 


L  - - - - 

Y  READ  INPUT  SPECIFICATIONS  OF  FILTER  TO  BE  DESIGNED 


a  READ (3, 6)  SPATE, FClOa,fCHU,H,kIPFLE,FSlOV,FSHIGH, STPLVl, IKINO, 
1ITtPE,uPOLy5,IPyOT,NPTS'ELOW,FhI5H 
f  FOrHAT ( 3F8.C,F6«0,eF 8, 0,F4,O, 211*212, Ia«2F; 3,0) 

IF ( ITTPE.Ea, 0)  GO  To  lOOO 


Y  aRITy  OUT  INPUT  SPECIFICATIONS  On  PRINTER 


»«lTE<»,7)  SRAfE , F yLOR. PCHI GH, m I PPlE , FsEOR,F SHI GH, STPLVL, I KINO, 
iitype*i,polls,ipeoi 

I  FOpRA T 1 1HI / 1 9X . ' If, PUT  sPECiFICATIONS'/V/lBX.  'SKATL  =  ',F8.2/V23X, 
,'FlLOV  s  ‘,F8.2//2uX,,fCHIuH  =  *  ,re.2//33X,  RIPPLE  =  • .F5.S//38X, 


Reproduced  from 
best  available  copy, 
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2'F«,L0*  =  '.Fa.2//4jX.'FSHI»H  s  ••Fa,2//‘teX.*STPLVL  =  '.F6.2//S3X. 
O'lKiNJ  =  '.U//58*,MTyPE  =  *.II//63x.'NP0lES  =  • » U//6SX*  •  IPLoT  = 

•*  sin 

FNYUsSrtATE/2. 

IF(FCU)X.(j).F.JY0.U«.FCnISH.ST.FNX8.0R.fSHIGH.6T.FNYU.0H,FSL0*.<,T. 
iFUYO)  GO  10  8 

IF<UPOuES.tO.O.ANu.STPL'<l..i.T.O.I  oO  To  8 
IF<Ik1uO.Ne.2.ANO.MPPlE.LI.0.|  Gy  TO  8 
1F(NPOlES.i<E.0.»NU.IKInU.Ew.3)  GO  TO  8 
IF ( ITYPE.to.ltANOsFCLOfc.GT ,FSLOX)  Go  TO  8 
lF(lTTHE.Eu.2.*NU.rCL0».LT.FSLO«)  GO  TO  8 

IF(lrrpE.t«.3.*ri0.(FCLy'Y.LI.FSo0X,0«.FcHlGH.GT.FSHlGH))  60  TO  8 
IF(lTrp£.Eu.4.AN0.lFCL0».Gr.FS4,0*,0K,F(.HIGH.LT,FSH!6H>)  60  TO  8 
00  TO  li 
o  «hiTE(4»9) 

f  FOKMAril8X,»0ESIG.<  IMPOSSIBLE-— -INCONSISTENT  SPECIFICATIONS') 
oC  TO  a 
Ii  lOkT inoe 


pU3,l‘tl59«:653!>8979324o0 

»CuT1=oTANiP1«FCL08/sRaTE> 

IF  ( ITTi'E.U)  ,3)  60  To  12 

»tur2=uJAN(pl«FCIIloH/S)<ATE) 

L>ir-Fs«CUT2-»CUTl 

Pi(oU=AtOll»»CUT2 

U  ir(liUi4U.Hu.2>  FPs8=(Io.0O0)»»(RIPPuE/10.0U0)-1.0U0 
IF  (UPJcES.NE.9.*iIi>.IKInO.Ns..3)  So  To  aO 


oETtuMlNE  M1.,IMU«  ORuER  FILTER  FhICM  MEETS  INPUT 
SPECIFICATIONS  UI.LESS  OHOER  IS  SPECIFIED 


*5 loPlaQTAui PI *F SLoW/SkATE  j 
uO  lo  (l3»i4.15.1a).lTYPE 

la  OHlGASoSTOpI/ACUTI 
00  To  18 

l o  o«u0A=»CUTi/XST0Pl 
00  To  IH 

la  »S)OP2=OTA„(Pl«FSMiGH/b!IATfc) 
lF<lTTp£.£o,4>  Go  TO  l6 

0Mt 4A=uA«a t ( USTOP l»»2-pK00 ) / ( «5T0P I «0 I FF ) > 
OML6AUOABa((9STOP2«*2-PPOy)/(»STOp2*OlFF) ) 
oo  To  17 

lo  0r,tuASoA6SU0STCPi«0lFF)/(P«00-xST0Pl**2)) 
OntOA  I  sOAOb  ( I  AST  OtV*OH  F I  /  ( PP00-»ST0P2»  »2 ) ) 

17  Im0NeOA1.lT. OMEGA)  O' EGASOPtGAl 

lo  00  TO  ( 19#s9*40) » iolND 

1 Y  YusSiHI ( (10,««(STPuVL/lO. )-l. I/EPS2) 

VU)  =  1.0 
V<i>sOwtGA 
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IF(OmEoa.<>|.V‘J>  G o  TO  2S 
DU  «!0  l=3««.l 

VU>s2,0«OwC6»*W|l“IJ-»II-2) 

IMV(I).LT.VN)  Gv  TO  *0 

NPoLES=I-l 

bo  TO  3? 

20  CONTINUE 

L 

21  «ikiTE(4.22( 

2c  F  okrtA  f  ( 40  A .  *  THE  SPuCIf  iCATlONS  REQUIRE  A  FILTtK  OE  ORDER  ,bT.  2P'» 
Gv  TO  3 
l 

23  kPoLESsl 
00  TO  30 

0 

SO  I.PuLESzALOolO  ( IP. **(STpLVL/10.)-l • )/( ?,*ALOOI OIOHEGA) )♦! 
IF(NpoLES.bT.20)  oO  To  21 
00  To  jO 
o 

<•0  AZi0.*»(STf'LVL/*0.  > 

BZuPSi/tA««-l.) 

C=1./(JKE0m«0.AECAI 
U3luLPiC(1.-U) 
tztLuPlu(i.-C) 

F=uLLPiC<dl 
OZlLlPTC (C ) 

UPuLESclF(r«E)/lO»o) 

IF  (UPOl.ES.  uT.  20 1  oO  To  21 
U=oE  ,<P  ( -NPoLES«P  I  »o/E) 

0-  .  • 

00  42  izl r iOO 
b=0Mll»l*i) 

IF (rt,Ll.l.u-30>  Go  TO  45 
Ozo+f' 

4*  COi.THhJE 
43  F=0. 

oo  47  1=1. iOO 
0=o**U»I) 

IF  (l  .L)  .I.u-30)  oo  TO  48 
K=f  *0 

47  COi.TINoE 

4o  o=(2.»o*»0,25«0/(1.+2.»F))«*4 
U=k  t-LPIC  ( l  ,-H) 

F=lLLP(C(B) 

L 

30  nzuPOLeS*  ( ( ITTPF  +11/2) 

»hiTEl4.51)  '1 

31  FohHAT(18*. 'THE  MiuIMUk  ORoE'I  FILTER  W..1CH  REE1S  THESE  SPEUFICATI 
.01.3  IS  ORUuR  '.lc. 

L 

L  CALLOLAIE  Thu  SEeOUU  O'IUORAN)  POlE  POSITIONS  IN  THE  S-PLANE 

L  OF  A  UNITY  okNURlOTM  LO*  PASS  FIlTER  OF  THE  REWIRED  ORDtR 

L  .(EAL12ING  in,.!  TmEIR  COMPLEX  CONJUGATES  ARE  ALSO  POLES 

V 

OO  I OuOzMoO (NpOLCS. 2 ) 

K=(IIP0lES*1)/2 
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fc*PNT=i.003/NPOLEa 
lF<lKI.,£).Ey.3)  GO  ro  90 

0«hGLE=PI*cXPNT 

A=1.0 

B=1.0 

IF ( 1k1nD.Eu.2)  Go  fO  ?o 

C=i<.  SwRT(l  300*1. 0u0/£pS2)*l. 0o0/UbaRT(EPS2.  )«*EXPNT 

A=.  500*  (C~l  .000/0 

B=.500«(C*1.000/C> 

70  ANgLE=0,0 

If  ( lOOu.EQ .  0 >  ANGlE=Da<''GLc./2  .  OD  0 
00  80  1=1.8 

pules ( i . l ) =-i . ooo* a«ocoS  < a„gle i 
POlES  ( 1 . 2 )  =B*DS  IN  (  angle  ) 

ANGLE=«NGLt*OANGLfc 
so  continue 

GO  TO  110 

90  *,i-(G/F)»uLOG(IObv.HT(l.+tl  '2> +1 , )/U$QkT (EPS2) ) 
fiusaHfd.-o 
INuIC=l 

-  92  (1=ELLPIC(F»F) 

&=tLLP|C(l.-F»F> 

G=utXP(-PUG/rf> 

ANgLE=PI»A/<2,»F) 

SN=0, 

00  94  JSO.J.COO 

0-(0**t  J*.3)  )/ll,-u»»<i;»J+l)  ) 

SUsSN*u*DSiN(  <2»J*1. )*ANGLt) 

IF(U.LI.l.u-JO)  GO  TO  9!, 

9*  CONTINUE 

9b  SN=(2.»PI«bN)/(f «o> 

GO  10  (96.48) .INDIl 
9u  »i<i=SN 

CNL=ObuRTU.-SNl  •  S.il) 

Di.l=DSuRT(i,-(F*SN).)***) 

IN0IC=2 

A=(l-IuDD)»E«EXPNf 
F  =oSQN |(C) 

1=0 

GO  10  92 
9o  SNo=SN 

CNO=DS  uRT  ( .  ,  •SU3*b.(3 ) 

ONj=DSgRT ( i . -C*SN J»SN3 ) 

1=1*1 

POlES  ( 1 . 1  >  =  < SNI  *C'U«CN J»ON 3 )  /  ( 1 .  -SI.I  »SN1  «PN3»UN  3 ) 
POlES ( I .2) =(SN3*DNi ) /( i.mSn1*SN1*0N3«0n3) 

IF (I, to, I. hND. 1000. EG, 1>  GO  TO  100 
2ExOS( J .11:0, 

ZEhOS  ( 1  •  2 )  51 .  /  (  OSG.tT  ( C )  *SN3 ) 

100  A=a*(2.*E*lXPNT> 

lFfl.LT. XI  GO  TO  92 


34 


TR  4629 


w  HAP  (HE  S-PuhNE  pOlES  TO  THfc  Z-Pl*NE  USING  TnE  BILINEAR 

t.  fRAHaFORMATloN  Ta*INo  THE  CUTOFF  FREQUENCY  AnO  TYPE  OF 

k.  ULTlP  INTO  mCCOuHT 

(.  ALSO,  IF  CLLlPTIo  IS  OESIHfO  CaR«T  OUT  THE  SAHE  procedure 

v  COR  1  HE  ZEKOs 

t 

110  1F<ITYP£.&(,2)  Go  TO  lTO 
J=1 

lFtlTYPE.U.I)  Gj  TO  i<*0 
IF  <  lOULi.ES. 0)  GO  |0  120 

Poe£Z <  i,l)  =  <POL£SU.l> +*10011  >/ (POleS (l,l)-*CUTt> 

Pol£Z( 1 ,2) -0,0 
lF(i.POLES.Ea.l)  uu  TO  300 
JS*. 

12U  L- j 
c 

UO  130  1SL,K 

Ai, -Ole^ ( 1 ,  D  +  jICUTI 

U=pOL£j(I,.:> 

UFpOLE  .>  ( 1 ,  i  )  -  ACUT 1 

POlLZ  ( u ,  1  )  =  (  A«0*U»i>  >  / 1  o*CHu«rt  > 

PulLZ  ( ii,  2 1  =  ( u»D-  A»U )  /  ( O»0*o»il  > 

Pul£Z(j*1.1>=P0L£Z<J,1) 

POlEZU*1>2)s-1,OoO*POleZ(  J,2) 

1F(IKI,,0.N1_.3)  GO  TO  125 

Ze^02( j, 1  )  =  (ZEROS!  1 , 2) +ZEKUS! 1 »2)“wCUTi*FCUT 1 )/ (Z£ROs! 1,2)* 
1ZEh05(1.2)*WCUT1«»CUT1> 

Z-.hOZ(  J,2)=2.  •ZEHua!l.  2)*WeUT1/ (ZEROS!  1>2>*Z£ROS(),Z)  +  i1CUT1»WCUT1> 
Zei,OZ(j+1,1>=ZEKOe(  J,I) 

ZEK0Z(O*l,«.)=-ZErtOAlJ,2) 

120  J=J»2 

uo  coi.tiool 

GO  TO  300 

L 

HO  C  =  1.0J0/«CoTl 

11 < 10UO.EU.0)  GO  ro  I*,0 
PolLZ ( 1 1 1 ) ilPOLES ( i , 1 ) »C )/ (C*P0LES( 1 , 1 ) > 

PoutZU* 21  =  0,0 
11  |HPOk.ES,^0,ll  Co  TO  300 
0=2 
ISO  L=o 

V 

Do  160  I =L , K 
A=hOLEs!1h)+C 
b=pOLEs<  I  ><. ) 

0=C-POuES(l,I) 

POi.F.Z!  j,l)i(A»0-8»o>/tu»D-K)*S> 

POLEZ(J,2)=(a*a*A»u)Z(0»D*o*R) 

PolEZI J*1<1)=P0LE4(J,1) 

POlEZI j*1,2)=-1.0o0*P0lLZ(J,Z) 

IF ( IkIuO.Ih.,3)  GO  to  155 

ZEh0Z(j,l)=(C»C-ZE«0S(I>21«ZeR0S(l,2))/(C*C+ZER0S(l,2>*ZEH05(I,2)> 
ZthOZ(j,2)=2.*ZFROJ(:,2)*C/(C*L*ZEHOS( 1,2) *ZEROS( 1,2)1 
ZEKOZ (o*l, 1) =Z£R0Z i J, 1 ) 

ZekOZ(j*1»£ )=-ZlRO*.( J,2) 

153  J=o*2 
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160  CONTINUE 
6v  TO  300 

170  1F(1TTP£.£«.3J  So  10  i90 
Ou  ISO  1=1. K 
A=POEEb( 1.1 ) 

POlES (  1 » 1 )=POLEF (1.1  )/(POLtSt  1 » 1 ! ••2+PoLES(  1,2 ) **2) 
POlES(1»2)=-POEESU.2)/(A«*2*POEE!»(1,2)»*2> 
1F(1x1nO.Nl,3)  60  TO  i<sO 

lF(l.EO.l.ANO.lODu.ca.l)  SO  TO  ISO 
A=ZLKOi,(  1.1) 

ZthO$(l.l)=ZEROS( 1.1)/ (ZEROS!  1 ,1  )**2+2t«0S(  1.2)  *»2) 
Zc.mCiS(1»2(=-ZERCS(..2)/(A**24ZER0S(1,2)**2) 

100  CONTINUE 

190  1NoIC=0 
190  CONTINUE 
J=1 

« 

DO  200  1=1, A, 

IT  I  1|Eu,2*nM}i1^Dui£0i  [)  u— U+2 

IP  l  (1.0T.2)  .OR.  <  1.4.0.2.AN0.10Uj,EO,0)  )  JSJ+4 

A=DlPF«DlFp»(P0Lt3ll,l)*»2-P0Lc.S(l,2)»».')-4.*PK0U 

b=i.»POLEStl.l)«PO4.ESU.2)»01FF«DlPF 

HMa&= ( A* A+o*B I **0.u9 

KAnSLcsDA  T  mNZ  ( b . A ) 

A=RMAo»DS1i«(RANGU./2.  ) 

8=nMA&»DC0 J ( A ANGLE/2 . ) 

C=i. FPA0D-U 1FF»P0LES (1,1) 

0=0 IFF  *P0Lc.S  (1.2) 

IF ( 1ND1C«E*.1)  60  TO  )9S 

INCUR 1=2 

IP { 1 .Eu. 1 *mND. IODu.Eo. | )  1NCRmT=1 
PulE2(u. 1 )=( ( 1.-PkuD*B>  «C-0«A ) / (C4C*0»0 ) 

PulEZ(»i.2)  =  )U*(1.-fR0D+i')4C*A)/(C»C+D»0) 
puuEz(J41n(.rrt.1)=(  (i.-prou-h)*c+o»A)/(C*c+o»u) 

POlEZ ( J* INlRNT , 2 ) s ( D« ( i » -PhOO-b ) *C*  A ) / ( C  *C  +U*D ) 
lPll.tu.l.AN&.lODu.EO.i)  oO  TO  200 
POlEZ(u*1»1)=°0*-EZ(J.1) 

PulEZ  ( w*  1  •  c )  =-POLc  a  ( J ,  2 ) 

POlEZ ( 0*3. i ) =POEEZ ( J+2 , 1 ) 

PulEZ ( J*3. * ! =-POLel  <  J+2 .2 ) 

60  TO  ZOO 

19o  1P(1,£w.1.mND.100u,EQ,i)  oO  To  200 

ZekOZ(J. 1)=( (l,-PAwn.B)»C-o*A>/(C*C+D40) 
ZtROZ(J.2)=(2»(l.-PROD*H)*L»A)/(C*C*D»U) 

ZEhOZ(  J+2.1)  =  (  (1. -PROD-0 )*t*D«A)/(0»C+Li»D) 
ZtROZ(u42.Z)=(D»(l.-PRoD-b)-C»A)/(C*C*0»0) 
ZEnOz(J+l.i)=ZE«OZ( 1.1) 

ZEmOZ( J*l«« )=-ZEROl(J.z) 

ZExOZ(u*3.1)=Z£ROZ(Jf2,1 ) 

ZEHOZ ( u+3>Z ) =-ZEROa (  J*2 . 2 ) 

ZOO  CONTINUE 

IFUKInO.Nl. 3. OP. 1NDIC»E8,11  SO  TO  300 

DO  210  1=1, K 
POl£S(1.1)=ZEROSU,1) 
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PUlES  <1 . 2  UZEROS  ( » .  2 ) 

210  continue 

llH/lC-l 
Go  To  14S 

30a  continue 

Mail 

1P(U1  (D.L„.J.ANJ.iODD.tu,0>  60  To  410 
IFiIkI..D.Eo.3)  NZl<ITyP£*1)/2 


o  P  ILL  UP  Tilt  lEROS  ARRAY  »ITH  THE  Z-PLANt  Z-lROg  (UNLESS  ELLIPTIC) 

t  IN  (TnlCH  CASl  THt  ZEROS  MATE  BEEN  CALCULATED. 


Do  '>00  1=1. NZ 

00  To  <310,32P.330.340).IT»PE 
3*0  ZlkOZU.Dl-I, 

ZEl.0ZW.2U0. 

Go  10  400 
)20  ZE„0ZU.IW1. 

ZerOZT  1.2W0. 

Go  10  hOO 

33U  Z«HOZa>l)  =  (-l.)»»(l+l) 

ZlrOZI I <2)=0. 
oO  TO  ->00 

340  Zti.02  <  1  >  1 )  =  1 1 .  -PRuu )  /  ( i .  4p«0ri  1 

ZL..OZU.2)-<-l.>*»(l  +  t>*2.»OSukTTPhOP>/(l.fPROU> 
“00  CONTINUE 


o  otTtnWINE  TIlL  COtFFKIENTS  IN  THt  DENOMINATOR  POLYNOMIAL 

o  OF  TnE  TRANSFER  fONCTIOn  BY  oOpBIMNG  THE  COMPLEX  CONJUGATE 

0  POLE  PAIRS  InTO  REAL  SECOND  order  FACTORS.  IF  THE  ORDER  IS 

v  oDO  IHEKt  WIlL  Ai.oO  oE  ONE  TERM  uF  ORDER  ONE.  CORRESPONDING 

o  TO  T.,£  HEAL  POLE 

0  FOLLo*  THt  SmME  PROCEDURE  TO  DETERMINE  THE  NUMERATOR  COEFFICIENTS 


410  u=* 

1F(I1YPE.G|.Z,oI>.1oDD.lP.O)  GO  TO  4S0 

CO|_t  ( i .  1  )=-POLEZ  ( i*  1  ( 

COlFU.2)='J. 

CulFNU.IU-ZEROZU.I) 

CulFn(i <2)=C. 

IF  (NPOLES.NE. I )  oO  To  440 
E=1.00U 

IF  ( 1  TY.-E  .Eo.  1 )  C=U . *C0EF  ( I .  I ) ) / ( 1 . FCotFNU  ,  1 ) ) 
IF  (  ITYpE.Eu.2)  C=U.-cOEF(1.1))/(1.-C0FFN<1.1)> 
GO  TO  oOO 
“40  U=._ 

“SO  L=w 


IF  1  ITYpE.gi  ,2)  Mr.iPOLtS 
Uo  4feO  10L.M 

COtF  U  .  1  i  =-P0LEZ  (J.l)  -,.0lEa  I  Jf  1 . 1 ) 

ColF ( 1 >2) =P0L£Z(U. 1 ) *PootZ( UFl . 1 ) FPOLEa ( U.2) •POLEZlU.2) 


37 


Tit  4629 


COtFN! i . 1 >  S-2EH0Z ( Ji 1 )-ZERoZ< J*1 • t ) 
COtFN<l#2>=2EROZ<J)l)*2EROX<J*l#l)+ZEROZlJ#2)*ZER02lJ*2> 

J=J*Z 

“60  CONTINUE 

L 

- - — - - - - - - - - ... - - - — 

C  HO*  (HAT  WE  NAVE  THE  POct-ZERO  PATTERN  OF  THE  DESIRED  TRANSFER 

C  FUNCTION)  HU)#  *E  HAVE  DETERMINED  HtZ)  UP  TO  SOME  CONSTANT. 

c  next  determine  That  constant. 

- - - - — - - - - ....... - .... 


Esi.OOO 

IKIOOu.EO.O.ANO.lKINO.NE.2)  E=1.0DO/OS8RTU  000+EPS21 
C=1.0D0 

00  TO  (470#480#490#470) # ITTPC 
c 

470  00  475  Izl-.M 

A?i#+CoEFNtI#l)+COcFN(I#2) 

6=1 . 000+C OuF  1 1 « 1 ) +COEF ( 1 . 2 ) 

C=c«<B/A) 

47a  Continue 
00  TO  500 
c 

460  DO  485  1=1. K 

A=1.-CoEFN(I.1)*COi_FN<I#2) 

B=1 , ODO-COtF ( I  #  1 )  K0EF  ( 1 , 2) 

C=l*(6/A) 

48o  CONTINUE 
GO  TO  500 
C 

490  THtTA=UATAN2 <2.»OSoRT!pHOO> . 1,-PRoU) 

F=dCGS(2.*)META) 

G=uC OS! THETA) 

0=USIN12.*|HETA> 

8=uSIN( THETA) 

DO  495  Isl,NPOLES 

B=uSQHr((FtC0EF(I.a)«G»C0EKI)2))»*2t(D*C0EF(I#l)*0)»*2) 
A=u5aHT((F+C0EFN(I,l)*6*C0tFN(  I  ,2))»»2+(0*C0EFN(I.1)*0)<>*2) 
C=c«(B/A) 

495  CONTINUE 
C 

500  CN<,T=E»C 
c 

C  CALCuLAT IONi  ARE  COMPLETE,  WRITE  OUT  THE  TRANSFER  FUNCTION 

C  AS  A  RATIO  OF  POlTNOmIALS  which  can  be  immediately 

C  IMPLEMENTED  iU  »C*SCACE'  form 

C 

WH1TE(4)510) 

510  format  (////ibx.'The  transfer  function, h|Z)#  of  the  oesired  filter 

i  Can  be  EXPRESSED' #/18x# 'IN  THE  FOLLOWING  FOR.-) 

If  ( lOOu.Efi.  I  .AND.  ITTPE.uT,3)  rRITc.(4#520) 

IF (IOOu.EO.O.OR.ITTPE.gT.2)  WRITE |4#530) 

520  FORMAT! //50X#‘ZFBU)')7X.'2w*2*B(2)»ZfBIJ)'*13X#'Z»*Z+B(N*a)*Z+B(N 

1!'./35X.'H(Z)  =  CNiT  *  -  •  - - -  *  ...  » 

2— /SoX>  UK!  1)'#7X,U*.2+C(2)»ZK!3>'.13X.'Z«»2+ 
4CCN-I)»ZKiN»'> 
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530  F0KMAT{//b0Xt'Z«*^*8(li»Z«ti<2)  '.7X,*Z»»2*B(3r  7»B(4)  •  •  13X. 

tN-jl»Z*0(N)»./3iX.'H(Zi  -  uNST  •  -  »  - - 

- -  •  ...  •  —,  r/50X,*Z«»2+C(l)»Z*C»2)'.7X 

j. *t**2*C (3)*Z»C (41 » . 13*. ‘2**2*C (N-i l*Z*C IN) ' ) 

•KlTEU.540)  CNST 

540  .FuhHAf  (//18X.»»HL.<E  Tilt  COEFFICIENTS  A»>E  GIVEN  BT',//3SX,'CNST  = 
i',u<3.13> 

INolC=U 

IKlOUu.EG.l  .ANU.lTTPl  .LT.31  IN0IC=1 
Liu  360  1  =  1, N 
Jl=(  I  +  illDIo  +  1 )  /2 

J2=2-MuO(  UINOIC  { l-l“OIC*  1 1 

»hiTE(4,55u)  I>Cuc.F(J1iJ2|  •  1  .CotFMol i J2> 

330  Fui.i'AI  (16a,  'C( '  Ilci  >  )  =  '>Ub3.15i10a>  'B(  '  i  12,  •  )=>, (1*3,15) 

360  toi.TIIIot 

6nl (£(4,5701 

670  FU„i  AT(///iPX,MHL  TRANSFER  FUNCTION  13  EXPRESSED  IN  TmE  AbCVE  FOR 
ih  so  1 MA T * / 18X ,  •  Tht  COtFFltlEllTS  CaN  Bl  USED  AS  GIVEN  IN  A  CASCADE 
A  inPLLi-CNTaTIO!.  OF  the  iILIER*) 

IF { lpLul.tu.O)  GO  TO  S 


C  .LOT  Tilt  NAGi.ITODt  AhO  PhASE  OF  H(Z1  AS  i  VARIES  On  1  HE 

i.  unit  circle  in  Tut  complex  plane 


»Lo*i2,*FL.a*PI/SHmTE 

SiU<iH=t.»F,.IGrt»PI/3R4Tt 

KAi,GE:„HlGr,-l>LOk 
UA= (F'HiGH-r  LOm) /1G ,0 
T  «oUXSt.»0A 
L=1 

DO  1.00  I  =  1,IOCO 

XL  L  1  =Ai_0*+i-L0AT  ( 1 )  *RANoE/I000 , 0 
SlMsSlMAU)  I 
S1i,2x  =  3IN(o.»X(I> ) 

COSXsCuS(XlD) 

Cosa'X=i.OS(a.»X<I>) 

W-=P. 

011,0 

IMlOOu.EG.O.OH.ITrHE.bT.Z)  GO  TO  3R0 
At  (COSA+CObFNL  1 , 1 )  I^.SINX*** 

U=tCOSx4COk.F  (1.  t>  >*«2*SlNX*«2 
0=uS9RUA/t) 

KIsColfMI.IUCOSa 
ht=COEF ( 1 , 1 1 +COSX 

if  t ahs lsina i ,lt , i ,u”io , and, ahs (Ri i,it,i,e-io)  ri=i.e-io 

lMAI|SLSINAl.Lr.l.i.-I0.ANO.Ar(S(R2)  .LT.1.E-I0)  H2=1.E-10 

KK=ATAi,2(SiNX,Rn-*TAN*(SINX.K*) 

Lt* 

IF  L'lROLtS.uO.l)  GO  (0  393 
360  1 01,11  Hot 

DO  590  J=L.M 

C=iCOSiX4CvEFU,  1)«COSa*COlFI  J,2>  l»*2*(SIN2X*CU-F( j»1)*SINX>**2 
AnCOSiXtCoEFiKJillACOSAFCuEFiHj,^) )  *•*♦  ISIN2afCOEFN(  J,  1)  aSIIIX  )  «*2 
B=g*03oRT(a/C) 

R1=S1i(EX4CoEFN(J.Ii»SIi,a 
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H*iC0S*X*Cu£FN(J.1)»C0s<*Co£FN( J.2) 

H3;S1i4.:X*Co£F(J.1)»S1Nx  v 
H4;COS2X*Cv£F(J,1)«COSx*COlF(J,2) 

lF|*aS(Rl).LT.l.E-i0,»NiJ.*DS(K2I.LT»l.t-10)  K2=l.t-10 
1F(AUL(R3).lT.!  .E-IO.An-I.AgSCHMI.lT.I.L-IO)  H4H.E-10 
NH3KR*aTAN,:(RI,'<2/-ATAn2<R3.R4) 

590  Cut.TiNuE 

*.>9s  i>- (U.li.o.ooi/chst)  so  ro  59o 
»(1>s2u.*DlOs10(CH;>T»B) 
lKmi.LT. -50.01  7(1)1-50.0 

uo  To  399 

59o  T ( 1 )3-d0, 

594  »(l)=X(l)»sR*Te/<2.0»Pi) 

NF4=lFlX(0.5»«N/Fll 

PiU(l):RR-t,»NPlAPi 

IF ( PHl (1)«lT,9.)  PHI ( i ) =Phl ( 1 )*2,aP1 
P,U  ( 11  =  180.  *PH1  ( 1 1/PJ 
hOe  COI.Tl.toE 

CA..L  MoDESuIZ.OI 

CAl*  OjJC Tu( 2, SOO.O.SOj. 0.3050 .0.3050,0) 

CAlL  5oSJtu(2,FLO«.-50.0,FrllS«,0.0) 

CalL  Gk1uG(2,DX.5.0,0.9,0.0) 

Cacl  LAbtn,(2,O.T«jL)X,O'0.2) 

CAlL  LaBELm(2, 1.10.0.0.5.1) 

C-lL  L1(4ESu(2.1000,X,7, 

CAli.  TULEuI  2.29,  'rREQuENC  1  II.  CTC.ES/T1ME  UNI  T  ' .».  'DEC  1BELS'  ,bt , 
I'PgOT  oF  T,.L  MAGn,  IUJE  OF  I  HE  FREuuENCT  RESPOHSl  OF  THt  FILTtR* ) 

CalL  OujC1u(2.30S0. .500. ,4095. ,  3050. ) 

Call  Sd9Jt„(2.0., 0. ,104b. .100.) 

Call  LlCNUu (2. 200*. 90. ,21, 'OE SCRIPT 1CN  OF  FIlTEHU 
GO  TO  (701.702.711.712|.ITTPE 
70i  CalL  Ll&NU„(2,2OO..0C..S.'lO«4  PASS'! 

GO  TO  / 0b 

70*  CAUL  Ll6NUw12.200..BO,,9,».iI6H  PASS') 

70s  CALL  LtGNU»(2.200.,60..e.'LUT0FF:  •) 

Call  HUMHKU (2,450. .60. , 7,2,FCL0*) 

IF  ( lKli.D.Eu.2)  Gu  TO  700 

CalL  LlGNO„ (2,200. , 55. 'H1PPl£=  »> 

CalL  NuMH(<j(2.450..55.  , a, 2, RIPPLE) 

Go  To  /SO 

711  call  LcGNOu (2, 200. , SO, ,0, 'oANUPASS' 1 
Go  TO  71b 

712  CAlL  Lc.GNOg  (2 , 200 ,  ,00, ,11,  *f)AN(,  RfoiCfl 
7)3  CalL  LlGNOu (2*200. ,&0,,V,'lUT0FF1s  •) 

CalL  LlGNOu( 2,200. , 55, , 4. 'CUT0FF2:  •) 

CALL  NoMBK j (2, 400. ,50,, 7, 2 ,FCL0X ) 

CALL  NoMBHv, ( 2 , 4 00  .  ,  S5 .  ,  7.2.FCH1SH) 

If  (1K1..U.E3.2!  GO  TO  750 

CAlL  Lw&NU^I 2*200, ,50. , B, 'a 1PPl£=  ') 

CAlL  NumBauI 2,450. ,50, ,4,2 .RIFPlE) 

7SU  GO  To  (001,802,003, ,1K|L0 

001  Call  LlC-HU«(2,200.,  75.  ,u,  'lMEotSHEv'  > 

GO  TO  050 

50*  CALL  LlGNUv,(2>  2C0.  ,7b, ,  1 1 ,  'BUT  T£7.«oRTH'  ) 

GO  TO  oSO 

803  CAlL  llGNU„(2,200..75.,o.'lLL1pT1C) 
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9S0  CALL  LtGNCy(Z,200.  ,70.  ,7,  *SRATE=  ') 

CAlL  nuM0Ro<Z»l*OO.»7o.,7,l,SRArE> 

C/.lL  LLGNO»(Z»2P0..65.,7,>uRt>£«=  •) 

CAi.L  NuM6Hv,(Z.»00.  ,65.  ,2,N) 

Call  PaGEUiz. 0.1.1, 

CAlL  OuJCtulZ, 500., 500., 4000,, 2600.) 

CAlL  SuejE»(Z,FLO».0.,FHIl,it,360,) 

CALL  u.<IDGiZ»DX#60..0.,0,| 
call  LaUELo<Z-0,T*u0X.0»8.2> 
call  LhSCLvIZ* t^^Oi f Oi^«i } 

CALL  L.NES»<Z.IOOO.X,P„I) 

CAlL  Ti;lE„<2«29, 'FREGulNCT  IN  CYClES/IIME  UNI T • .7, 'DEGREES’ >57. 'P 
1LOT  OF  IKE  PHASE  ur  THL  FRlQUEnCY  RESPONSE  OF  Tut  FILTER') 

Call  PmGEGiZ, 0.1.1, 


L  lOMPuTE  AW)  PLOT  Tht  UNIT  SAMPLE  RESPONSt  AS  THE  INVERSE 

L  Z-TRaNSFOHN  uF  H(Z) 

C  ------ — - - 

I 

Uu  9o0  1=1. N 
K*=POLlZ(1,I> 

RZ=P0L«.Z<I,21 

,<3=ZERvZ<I.ll 

rt4iZlrtuZ<I.2, 

ZPuLES 1 I ) =tMPL  X ( HI . H2 ) 

ZTt''P<il=l. 

<00  ZZl'<05(  II=tMPLX<R.).H4) 


<U0 

920 


UU  920  1=1, N 

PTlMPIlCMPvX<1..0.) 

P1l”P«!=CMP1.XU.»9.) 

UU  910  J=t,N 

PIlMP1-:PTEmP1«<7PulES<  1 1-Z2ERUSU) ) 
IHJ.tu.D  GO  To  vlo 
P I L  NP2  =PTt ,.,P2  *  ( ZPuvES  1 1 )  -ZPOLlS  ( J ) ) 
CUi.TlNuL 

UL)AU|=PT,.MP1/PTli-.P2 

CUi.TINuE 


Hi t IT-O, 

H2ll,-'  NST 
HuaX^vijsT 

inipis.ta.c)  npij=ioo 

P!  5=FLUAT  IlPT.') 

UO  9b0  1=2, NP.‘  > 

ZTL«Pi=0. 
uu  9jn  jsi.m 

ZlLMPi:2TtMP2*rtrTrt(J)«2lEMPlJI 

2ItAP(J)=Z)EHP(J)»AP0LtitJ> 

930  CO».T  INut 

HI  ( 1 , -FLOA ( ( I— I ) 

(«(l)=vNST.RCALtZlL9P2| 

IF  t  AHb(H2(l>  i  .GT.H.iAX)  HMAX=AuS(M2U  )) 
9b0  CONTINUE 

HF.AA=J.05*rlM*X 
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II-  (NMAn.LT. 0.002) 
HMlUA-NMAX 
DllHMAx/2.;, 
OASI)PTi/10 


HMAa-0.002 


CalL  Ou-JC  Ti-.Z,  BOO..  SCO.! ‘4000.  i,;800.) 

CAtL  '.je.‘.«(2»0.>ii:4lM,(,IS.rt*)AA) 

CmlL  0  "JGU'PTS'HhAX'U.'n,  ) 

CmlL  l«BLLg<Z.O>Im,0.4) 

CAlL  LmB£Lg(Z# 1 # Of  #0.6. 3) 

CALL  L1N£S«(2.NPTs.H1,h2) 

CAlL  TlTLEwU.lS.  'aESPoNS£  NUMbER*  *12»  *  OUTPUT  VAlO£’  .38.  'PLOT  OF  U 
iNIT  SAnPLE  RESPONSi  OF  HLTEB'l 
call  PaGEG(Z,0<1i1| 

GO  10  3 


4000  COi.IIUmE 

CAlL  EaITGiZ) 
tl<L 


LflU  CF  LCMi  lLATl jlC 
SlocS  SfMoOLlC. 

S1od9  COul  AEluCAjAhLl 


r.O  DIAGNOSTICS. 


23 

JUL 

73 

08:<ti:i4 

0 

00246670 

14 

711 

23 

JUL 

73 

oa:4i:i4 

l 

03272232 

46 

1 

0 

03272312 

14 

347 

(DELETED) 

(OELE’EO) 


G*(.  AGl  iioby 


Ia:481J«.  SO 
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